function [Outp,Exit,Addf,Dcy] = simulate(varargin)
% simulate  Simulate model.
%
%
% Syntax
% =======
%
%     S = simulate(M,D,Range,...)
%     [S,Flag,AddF,Delta] = simulate(M,D,Range,...)
%
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object.
%
% * `D` [ struct | cell ] - Input database or datapack from which the
% initial conditions and shocks from within the simulation range will be
% read.
%
% * `Range` [ numeric | char ] - Simulation range.
%
%
% Output arguments
% =================
%
% * `S` [ struct | cell ] - Database with simulation results.
%
%
% Output arguments in nonlinear simulations
% ===========================================
%
% * `ExitFlag` [ cell | empty ] - Cell array with exit flags for
% nonlinearised simulations.
%
% * `AddF` [ cell | empty ] - Cell array of tseries with final add-factors
% added to first-order approximate equations to make nonlinear equations
% hold.
%
% * `Delta` [ cell | empty ] - Cell array of tseries with final
% discrepancies between LHS and RHS in equations marked for nonlinear
% simulations by a double-equal sign.
%
%
% Options
% ========
%
% * `'anticipate='` [ *`true`* | `false` ] - If `true`, real future shocks
% are anticipated, imaginary are unanticipated; vice versa if `false`.
%
% * `'contributions='` [ `true` | *`false`* ] - Decompose the simulated
% paths into contributions of individual shocks.
%
% * `'dbOverlay='` [ `true` | *`false`* | struct ] - Use the function
% `dboverlay` to combine the simulated output data with the input database,
% (or a user-supplied database); both the data preceeding the simulation
% range and after the simulation range are appended.
%
% * `'deviation='` [ `true` | *`false`* ] - Treat input and output data as
% deviations from balanced-growth path.
%
% * `'dTrends='` [ *`@auto`* | `true` | `false` ] - Add deterministic
% trends to measurement variables.
%
% * `'ignoreShocks='` [ `true` | *`false`* ] - Read only initial conditions
% from input data, and ignore any shocks within the simulation range.
%
% * `'method='` [ *`'firstorder'`* | `'selective'` | `'global'` ] - Method
% of running simulations; `'firstorder'` means first-order approximate
% solution (calculated around steady state); `'selective'` means
% equation-selective nonlinear method; `'global'` means global nonlinear
% method (available only in models with no leads).
%
% * `'plan='` [ plan ] - Specify a simulation plan to swap endogeneity and
% exogeneity of some variables and shocks temporarily, and/or to simulate
% some nonlinear equations.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the
% command window.
%
% * `'sparseShocks='` [ `true` | *`false`* ] - Store anticipated shocks
% (including endogenized anticipated shocks) in sparse array.
%
%
% Options for equation-selective nonlinear simulations
% =====================================================
%
% * `'solver='` [ *`@qad`* | `@fsolve` | `@lsqnonlin` ] - Solution
% algorithm; see Description.
%
% * `'maxNumelJv='` [ numeric | *`1e6`* ] - Maximum number of data points
% (nonlinear plus exogenized) allowed for a nonrecursive algorithm in the
% nonlinear equation updating step; if exceeded, a recursive
% (period-by-period) simulation is used to update nonlinear equations
% instead.
%
% * `'nonlinPer='` [ numeric | *`@all`* ] - Horizon (number of periods from
% the beginning of the simulation, and from the beginning of each
% simulation segment) over which nonlinearities will be preserved; the
% remaining periods will be simulated using first-order approximate
% solution.
%
%
% Options for equation-selective nonlinear simulations with @qad solver
% =======================================================================
%
% * `'addSstate='` [ *`true`* | `false` ] - Add steady state levels to
% simulated paths before evaluating nonlinear equations; this option is
% used only if `'deviation=' true`.
%
% * `'display='` [ *`true`* | `false` | numeric | Inf ] - Report iterations
% on the screen; if `'display=' N`, report every `N` iterations; if
% `'display=' Inf`, report only final iteration.
%
% * `'error='` [ `true` | *`false`* ] - Throw an error whenever a
% nonlinear simulation fails converge; if `false`, only an warning will
% display.
%
% * `'lambda='` [ numeric | *`1`* ] - Initial step size (between `0` and
% `1`) for add factors added to nonlinearised equations in every iteration;
% see also `'nOptimLambda='`.
%
% * `'nOptimLambda='` [ numeric | `false` | *`1`* ] - Find the optimal step
% size on a grid of 10 points between 0 and `'lambda='` before each of the
% first `'nOptimLambda='` iterations; if `false`, the value assigned to
% `Lambda` is used and no grid search is performed.
%
% * `'reduceLambda='` [ numeric | *`0.5`* ] - Reduction factor (between `0`
% and `1`) by which `lambda` will be multiplied if the nonlinear
% simulation gets on an divergence path.
%
% * `'upperBound='` [ numeric | *`1.5`* ] - Multiple of all-iteration
% minimum achieved that triggers a reversion to that iteration and a
% reduciton in `lambda`.
%
% * `'maxIter='` [ numeric | *`100`* ] - Maximum number of iterations.
%
% * `'tolerance='` [ numeric | *`1e-5`* ] - Convergence tolerance.
%
%
% Options for nonlinear simulations with Optim Tbx solver
% ========================================================
%
% * `'optimSet='` [ cell | struct ] - Optimization Tbx options.
%
%
% Options for global nonlinear simulations
% =========================================
%
% * `'optimSet='` [ cell | struct ] - Optimization Tbx options.
%
% * `'solver='` [ `@fsolve` | *`@lsqnonlin`* ] - Solution algorithm; see
% Description.
%
%
% Description
% ============
%
% The function `simulate(...)` simulates a model on the specified
% simulation range. By default, the simulation is based on a first-order
% approximate solution (calculated around steady state). To run nonlinear
% simulations, use the option `'nonlinear='` (to set the number of periods
%
% Output range
% -------------
%
% Time series in the output database, `S`, are are defined on the
% simulation range, `Range`, plus include all necessary initial conditions,
% ie. lags of variables that occur in the model code. You can use the
% option `'dboverlay='` to combine the output database with the input
% database (ie. to include a longer history of data in the simulated
% series).
%
%
% Deviations from steady-state and deterministic trends
% ------------------------------------------------------
%
% By default, both the input database, `D`, and the output database, `S`,
% are in full levels and the simulated paths for measurement variables
% include the effect of deterministic trends, including possibly exogenous
% variables. The default behavior can be changed by changing the options
% `'deviation='` and `'dTrends='`.
%
% The default value for `'deviation='` is false. If set to `true`, then the
% input database is expected to contain data in the form of deviations from
% their steady state levels or paths. For ordinary variables (ie. variables
% whose log status is `false`), it is $x_t-\Bar x_t$, meaning that a 0
% indicates that the variable is at its steady state and e.g. 2 indicates
% the variables exceeds its steady state by 2. For log variables (ie.
% variables whose log status is `true`), it is $x_t/\Bar x_t$, meaning that
% a 1 indicates that the variable is at its steady state and e.g. 1.05
% indicates that the variable is 5 per cent above its steady state.
%
% The default value for `'dTrends='` is `@auto`. This means that its
% behavior depends on the option `'deviation='`. If `'deviation=' false`
% then deterministic trends are added to measurement variables, unless you
% manually override this behavior by setting `'dTrends=' false`.  On the
% other hand, if `'deviation=' true` then deterministic trends are not
% added to measurement variables, unless you manually override this
% behavior by setting `'dTrends=' true`.
%
%
% Simulating contributions of shocks
% -----------------------------------
%
% Use the option `'contributions=' true` to request the contributions of
% shocks to the simulated path for each variable; this option cannot be
% used in models with multiple alternative parameterizations or with
% multiple input data sets.
%
% The output database, `S`, contains Ne+2 columns for each variable, where
% Ne is the number of shocks in the model:
%
% * the first columns 1...Ne are the
% contributions of the Ne individual shocks to the respective variable;
%
% * column Ne+1 is the contribution of initial condition, th econstant, and
% deterministic trends, including possibly exogenous variables;
%
% * column Ne+2 is the contribution of nonlinearities in nonlinear
% simulations (it is always zero otherwise).
%
% The contributions are additive for ordinary variables (ie. variables
% whose log status is `false`), and multplicative for log variables (ie.
% variables whose log status is `true`). In other words, if `S` is the
% output database from a simulation with `'contributions=' true`, `X` is an
% ordinary variable, and `Z` is a log variable, then
%
%     sum(S.X,2)
%
% (ie. the sum of all Ne+2 contributions in each period, ie. summation goes
% across 2nd dimension) reproduces the final simulated path for the
% variable `X`, whereas
%
%     prod(S.Z,2)
%
% (ie. the product of all Ne+2 contributions) reproduces the final
% simulated path for the variable `Z`.
%
%
% Simulations with multiple parameterisations and/or multiple data sets
% ----------------------------------------------------------------------
%
% If you simulate a model with `N` parameterisations and the input database
% contains `K` data sets (ie. each variable is a time series with `K`
% columns), then the following happens:
%
% * The model will be simulated a total of `P = max(N,K)` number of times.
% This means that each variables in the output database will have `P`
% columns.
%
% * The 1st parameterisation will be simulated using the 1st data set, the
% 2nd parameterisation will be simulated using the 2nd data set, etc. until
% you reach either the last parameterisation or the last data set, ie.
% `min(N,K)`. From that point on, the last parameterisation or the last
% data set will be simply repeated (re-used) in the remaining simulations.
%
% * Formally, the `I`-th column in the output database, where `I = 1, ...,
% P`, is a simulation of the `min(I,N)`-th model parameterisation using the
% `min(I,K)`-th input data set number.
%
%
% Equation-selective nonlinear simulations
% -----------------------------------------
%
% The equation-selective nonlinear simulation approach is invoked by
% setting `'method=' 'selective'`. In equation-selective nonlinear
% simulations, the solver tries to find add-factors to user-selected
% nonlinear equations (ie. equations with `=#` instead of the equal sign in
% the model file) in the first-order solution such that the original
% nonlinear equations hold for simulated trajectories (with expectations
% replaced with actual leads).
%
% Two numerical approaches are available, controlled by the option
% `'solver='`:
%
% * '`QaD`' - a quick-and-dirty, but less robust method (default);
%
% * `@fsolve`, `@lsqnonlin` - which are standard Optimization Tbx routines,
% slower but likely to converge for a wider variety of simulations.
%
%
% Global nonlinear simulations
% -----------------------------
%
% The global nonlinear simulation approach is invoked by setting `'method='
% 'global'` and is available only in models with no leads (expectations).
% In global nonlinear simulations, the entire model is solved as a system
% of nonlinear equations, period by period, using one of the following two
% Optimization Tbx routines: `@fsolve` or `@lsqnonlin` (default).
%
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

[This,Inp,Range,varargin] = ...
    irisinp.parser.parse('model.simulate',varargin{:});
opt = passvalopt('model.simulate',varargin{:});

if ischar(opt.solver) && strcmpi(opt.solver,'plain')
    opt.solver = @qad;
end

%--------------------------------------------------------------------------

% Input struct to the backend functions in `+simulate` package.
s = struct();

ny = sum(This.nametype == 1);
nx = size(This.solution{1},1);
nb = size(This.solution{1},2);
nf = nx - nb;
ne = sum(This.nametype == 3);
ng = sum(This.nametype == 5);
nAlt = size(This.Assign,3);
[~,maxLead] = mymaxshift(This);
isBkw = maxLead == 0;

nPer = length(Range);
s.NPer = nPer;

% Simulation plan.
isSwap = isa(opt.plan,'plan') ...
    && nnzendog(opt.plan) > 0 && nnzexog(opt.plan) > 0;

% Get initial condition for alpha.
% alpha is always expanded to match nAlt within datarequest(...).
[alp0,x0,nanInit] = datarequest('init',This,Inp,Range);
if ~isempty(nanInit)
    if isnan(opt.missing)
        nanInit = unique(nanInit);
        utils.error('model:simulate', ...
            'This initial condition is not available: ''%s''.', ...
            nanInit{:});
    else
        alp0(isnan(alp0)) = opt.missing;
    end
end
nInit = size(alp0,3);

% Get shocks; both reals and imags are checked for NaNs within
% datarequest(...).
if ~opt.ignoreshocks
    eInp = datarequest('e',This,Inp,Range);
    % Find the last anticipated shock to determine t+k for expansion.
    if opt.anticipate
        lastEa = utils.findlast(real(eInp));
    else
        lastEa = utils.findlast(imag(eInp));
    end
    nShock = size(eInp,3);
else
    lastEa = 0;
    nShock = 0;
end
s.LastEa = lastEa;

% Check for option conflicts.
doChkConflicts();

yTune = [ ];
xTune = [ ];
lastEndgA = 0;
lastEndgU = 0;
nSwap = 0;
if isSwap
    % Get anchors; simulation range and plan range must be identical.
    % Get data for measurement and transition tunes.
    doGetPlanData();
else
    s.Anch = false(ny+nx+ne+ne,nPer);
    s.Wght = sparse(ne+ne,nPer);
end
s.LastEndgA = lastEndgA;
s.LastEndgU = lastEndgU;

% Get exogenous variables in dtrend equations.
G = datarequest('g',This,Inp,Range);
nExog = size(G,3);

% Total number of cycles.
nLoop = max([1,nAlt,nInit,nShock,nSwap,nExog]);
s.NLoop = nLoop;

Exit = cell(1,nLoop);
Addf = cell(1,nLoop);
Dcy = cell(1,nLoop);

s = myprepsimulate1(This,s,opt);
doChkNonlinConflicts();

% Initialise handle to output data.
xRange = Range(1)-1 : Range(end);
if ~opt.contributions
    hData = hdataobj(This,xRange,nLoop);
else
    hData = hdataobj(This,xRange,ne+2,'Contributions=',@shock);
end

% Preallocate array for time-varying parameter updates.
if s.IsUpdate
    nUpd = length(s.Update.PosUpdate);
    pData = nan(nUpd,nPer,nLoop);    
end

% Main loop
%-----------
if opt.progress ...
        && ( This.IsLinear ...
        || isequal(opt.display,0) || isequal(opt.display,'off') ...
        || isequal(opt.display,'none') )
    s.progress = progressbar('IRIS model.simulate progress');
else
    s.progress = [];
end

for iLoop = 1 : nLoop
    s.ILoop = iLoop;
    
    if iLoop <= nAlt
        % Update solution and other loop-dependent info to be used in this
        % simulation round.
        s = myprepsimulate2(This,s,iLoop);
    end

    % Get current initial condition for the transformed state vector,
    % current shocks, and measurement and transition tunes.
    doGetData();
    
    % Compute deterministic trends if requested. Do not compute the dtrends
    % in the `+simulate` package because they are dealt with differently when
    % called from within the Kalman filter.
    dTnd = [ ];
    if ny > 0 && opt.dtrends
        dTnd = mydtrendsrequest(This,'range',Range,s.G,iLoop);
        if isSwap
            % Subtract deterministic trends from measurement tunes.
            s.Tune(1:ny,:) = s.Tune(1:ny,:) - dTnd;
        end
    end
    
    % Call the backend package `simulate`
    %-------------------------------------
    exit = [ ];
    dcy = [ ];
    addf = [ ];
    s.y = [ ]; % Measurement variables.
    s.w = [ ]; % Transformed transition variables, w := [xf;alpha].
    s.v = [ ]; % Correction vector for nonlinear equations.
    s.M = [ ];
    
    if s.IsUpdate
        % Simulate period by period, update steady state in each period.
        [s.y,s.xx,s.Ea,s.Eu,s.p] = mysimulateper(This,s);
        pData(:,:,iLoop) = s.p;
    else
        switch s.Method
            case 'firstorder'
                % Linear simulations.
                if opt.contributions
                    if isSwap
                        [s.y,s.xx,s.Ea,s.Eu] = simulate.linear.run(s,Inf);
                    end
                    [s.y,s.xx,s.Ea,s.Eu] = ...
                        simulate.linear.contributions(s,Inf);
                else
                    [s.y,s.xx,s.Ea,s.Eu] = simulate.linear.run(s,Inf);
                end
            case 'selective'
                % Equation-selective nonlinear simulations.
                if opt.contributions
                    % Simulate linear contributions of shocks.
                    c = struct();
                    [c.y,c.xx,c.Ea,c.Eu] = ...
                        simulate.linear.contributions(s,Inf);
                end
                % Simulate contributions of nonlinearities residually.
                [s.y,s.xx,s.Ea,s.Eu,~,exit,dcy,addf] = ...
                    simulate.selective.run(s);
                if opt.contributions
                    c.y(:,:,ne+2) = s.y - sum(c.y,3);
                    c.xx(:,:,ne+2) = s.xx - sum(c.xx,3);
                    s.y = c.y;
                    s.xx = c.xx;
                    s.Ea = c.Ea;
                    s.Eu = c.Eu;
                end
            case 'global'
                % Simulate period by period.
                [s.y,s.xx,s.Ea,s.Eu] = mysimulateper(This,s);
        end
    end
    %-------------------------------------------------------------------
    % Beyond this point, only `s.y`, `s.xx`, `s.Ea` and `s.Eu` are used
    
    % Diagnostics output arguments for selective nonlinear simulations.
    if isequal(s.Method,'selective');
        Exit{iLoop} = exit;
        Dcy{iLoop} = dcy;
        Addf{iLoop} = addf;
    end
    
    % Add measurement detereministic trends.
    if ny > 0 && opt.dtrends
        % Add to trends to the current simulation; when `'contributions=' true`, we
        % need to add the trends to (ne+1)-th simulation (ie. the contribution of
        % init cond and constant).
        if opt.contributions
            s.y(:,:,ne+1) = s.y(:,:,ne+1) + dTnd;
        else
            s.y = s.y + dTnd;
        end            
    end

    % Assign output data.
    doAssignOutp();
    
    % Add equation labels to add-factor and discrepancy series.
    if isequal(s.Method,'selective') && nargout>2
        label = s.Selective.EqtnLabelN;
        Dcy{iLoop} = permute(Dcy{iLoop},[2,1,3]);
        Dcy{iLoop} = tseries( Range(1), Dcy{iLoop}, label );
        Addf{iLoop} = permute(Addf{iLoop},[2,1,3]);
        nSgm = size(Addf{iLoop},3);
        label = repmat(label,1,1,nSgm);
        Addf{iLoop} = tseries( Range(1), Addf{iLoop}, label );
    end

    % Update progress bar.
    if ~isempty(s.progress)
        update(s.progress,s.ILoop/s.NLoop);
    end
end % for

% Post mortem
%-------------
if isSwap
    % Throw a warning if the system is not exactly determined.
    doChkDetermined();
end

% Convert hdataobj to struct. The comments assigned to the output series
% depend on whether contributions=true or false.
Outp = hdata2tseries(hData);
if s.IsUpdate
    doParamUpdDb();
end

% Overlay the input (or user-supplied) database with the simulation
% database.
if isequal(opt.dboverlay,true)
    Outp = dboverlay(Inp,Outp);
elseif isstruct(opt.dboverlay)
    Outp = dboverlay(opt.dboverlay,Outp);
end


% Nested functions...


%**************************************************************************


    function doChkNanExog()
        % Check for NaNs in exogenised variables.
        inx1 = s.Anch(1:ny+nx,:);
        inx2 = [any(~isfinite(yTune),3);any(~isfinite(xTune),3)];
        inx3 = [any(imag(yTune) ~= 0,3);any(imag(xTune) ~= 0,3)];
        inx = any(inx1 & (inx2 | inx3),2);
        if any(inx)
            list = myvector(This,'yx');
            utils.error('model:simulate', ...
                ['This variable is exogenised to NaN, Inf or ', ...
                'complex number: ''%s''.'], ...
                list{inx});
        end
    end % doChkNanExog()


%**************************************************************************


    function doChkDetermined()
        if nnzexog(opt.plan) ~= nnzendog(opt.plan)
            utils.warning('model:simulate', ...
                ['The number of exogenised data points (%g) does not ', ...
                'match the number of endogenised data points (%g).'], ...
                nnzexog(opt.plan),nnzendog(opt.plan));
        end
    end % doChkDetermined()


%**************************************************************************


    function doAssignOutp()
        n = size(s.xx,3);
        % Add pre-sample init cond to x.
        xf = [nan(nf,1,n),s.xx(1:nf,:,:)];
        xb = s.xx(nf+1:end,:,:);
        if opt.contributions
            pos = 1 : ne+2;
            xb = [zeros(nb,1,ne+2),xb];
            xb(:,1,ne+1) = x0(:,1,min(iLoop,end));
            g = zeros(ng,nPer,ne+2);
            g(:,:,ne+1) = s.G;
        else
            pos = iLoop;
            xb = [x0(:,1,min(iLoop,end)),xb];
            g = s.G;
        end
        % Add current results to output data.
        if opt.anticipate
            e = s.Ea + 1i*s.Eu;
        else
            e = s.Eu + 1i*s.Ea;
        end
        hdataassign(hData,pos, ...
            { ...
            [nan(ny,1,n),s.y], ...
            [xf;xb], ...
            [nan(ne,1,n),e], ...
            [], ...
            [nan(ng,1,n),g], ...
            });
    end % doAssignOutput()


%**************************************************************************


    function doChkConflicts()
        % The option 'contributions=' option cannot be used with the 'plan='
        % option, with multiple parameterisations, or multiple data sets.
        if opt.contributions
            if nAlt > 1 || nInit > 1 || nShock > 1
                utils.error('model:simulate', ...
                    ['Cannot simulate(...) ', ...
                    'models with multiple alternative parameterizations ', ...
                    'with option contributions=true.']);
            end
            if nInit > 1 || nShock > 1
                utils.error('model:simulate', ...
                    ['Cannot simulate(...) ', ...
                    'multiple data sets ', ...
                    'with option contributions=true.']);
            end
        end
    end % doChkConflicts()


%**************************************************************************


    function doChkNonlinConflicts()
        if isequal(s.Method,'selective') ...
                && lastEndgU > 0 && lastEndgA > 0
            utils.error('model:simulate', ...
                ['Cannot simulate(...) with option method=selective and ', ...
                'both anticipated and unanticipated exogenized shocks.']);
        end
        if isequal(s.Method,'selective') && s.IsUpdate
            utils.error('model:simulate', ...
                ['Cannot simulate(...) with option method=selective and ', ...
                'sstateUpdate=true.']);
        end
        if s.IsUpdate && s.IsDeviation
            utils.error('model:simulate', ...
                ['Cannot simulate(...) ', ...
                'with options deviation=true ', ...
                'and sstateUpdate=true.']);
        end
        if ~isBkw && isequal(s.Method,'global')
            utils.error('model:simulate', ...
                ['Cannot simulate(...) forward-looking models ', ...
                'with option method=global.']);
        end
    end % doChkNonlinConflicts()


%**************************************************************************


    function doGetData()        
        % Get current initial condition for the transformed state vector,
        % and current shocks.
        s.x0 = x0(:,1,min(iLoop,end));
        s.Alp0 = alp0(:,1,min(iLoop,end));
        if opt.ignoreshocks
            s.Ea = zeros(ne,nPer);
            s.Eu = zeros(ne,nPer);
        else
            if opt.anticipate
                s.Ea = real(eInp(:,:,min(iLoop,end)));
                s.Eu = imag(eInp(:,:,min(iLoop,end)));
            else
                s.Ea = imag(eInp(:,:,min(iLoop,end)));
                s.Eu = real(eInp(:,:,min(iLoop,end)));
            end
        end
        if opt.sparseshocks
            s.Ea = sparse(s.Ea);
        end
        % Current tunes on measurement and transition variables.
        if isSwap
            s.Tune = [ ...
                yTune(:,:,min(iLoop,end)); ...
                xTune(:,:,min(iLoop,end)); ...
                ];
        else
            s.Tune = sparse(ny+nx,nPer);
        end
        % Exogenous variables in dtrend equations.
        s.G = G(:,:,min(iLoop,end));
    end % doGetData()


%**************************************************************************


    function doGetPlanData()
        [yAnch,xAnch,eaAnch,euAnch,~,~,eaWght,euWght] = ...
            myanchors(This,opt.plan,Range,opt.anticipate);
        s.Anch = [yAnch;xAnch;eaAnch;euAnch];
        s.Wght = [eaWght;euWght];
        % Get values for exogenised data points.
        if any(yAnch(:))
            yTune = datarequest('y',This,Inp,Range);
        else
            yTune = zeros(ny,nPer);
        end
        if any(xAnch(:))
            xTune = datarequest('x',This,Inp,Range);
        else
            xTune = zeros(nx,nPer);
        end
        % Check for NaNs in exogenised variables.
        doChkNanExog();
        nSwap = max(size(yTune,3),size(xTune,3));
        lastEndgA = utils.findlast(eaAnch);
        lastEndgU = utils.findlast(euAnch);
    end % doGetSwapData()


%**************************************************************************


    function doParamUpdDb()
        temp = tseries();
        for ii = 1 : length(s.Update.PosUpdate)
            iiName = This.name{s.Update.PosUpdate(ii)};
            iiData = permute(pData(ii,:,:),[2,3,1]);
            Outp.(iiName) = replace(temp,iiData,range(1));
        end
    end % doParamUpdDb


end % main
